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An analog to the equations of compressible flow that is based on the inviscid Burgers equation 
is utilized to investigate the effect of spatial discreteness of energy release on the propagation of 
a detonation wave. While the traditional Chapman-Jouguet (CJ) treatment of a detonation wave 
assumes that the energy release of the medium is homogeneous through space, the system examined 
here consists of sources represented by 5-functions embedded in an otherwise inert medium. The 
sources are triggered by the passage of the leading shock wave following a delay that is either 
of fixed period or randomly generated. The solution for wave propagation through a large array 
(10^-10^) of sources in one dimension can be constructed without the use of a finite difference 
approximation by tracking the interaction of sawtooth-profiled waves for which an analytic solution 
is available. A detonation-like wave results from the interaction of the shock and rarefaction waves 
generated by the sources. The measurement of the average velocity of the leading shock front for 
systems of both regular, fixed-period and randomized sources is found to be in close agreement with 
the velocity of the equivalent CJ detonation in a uniform medium wherein the sources have been 
spatially homogenized. This result may have implications for the applicability of the CJ criterion 
to detonations in highly heterogeneous media (e.g., polycrystalline, solid explosives) and unstable 
detonations with a transient and multidimensional structure (e.g., gaseous detonation waves). 


I. INTRODUCTION 

Detonation waves are supersonic combustion waves in 
which a shock, propagating through a reactive medium, 
initiates energy release in the medium that, in turn, 
supports the initial shock wave. The coupled shock and 
reaction complex is nearly always observed to propagate 
at a velocity predicted by the Chapman-Jouguet (CJ) 
criterion, which states that the flow, upon reaching the 
end of the reaction zone, moves away from the wave at 
the local sonic velocity (i.e., Mach number unity) in the 
wave-fixed reference frame. This criterion is equivalent 
to stating that the leading edge of the rarefaction 
wave in the unsteady expansion of detonation products 
downstream of the wave is equal to the velocity of the 
detonation wave itself. The CJ criterion provides the 
necessary condition to close the conservation laws of 
mass, momentum, and energy and thus solve for the det¬ 
onation velocity and thermodynamic state of the reacted 
medium, provided the equation of state is known. The 
success of this simple criterion, which was formulated 
more than a century ago mm, is intriguing since it is 
formulated for a wave propagating at a steady velocity 
into a uniform medium, a scenario that is almost never 
encountered in experiments examining detonation waves. 
In gaseous explosives, the structure of the detonation 
front is always unstable and multidimensional, consisting 
of detonation cells defined by transverse shock waves 
that propagate across the front. [5] Similar features are 
believed to exist in liquid explosives as well.|l] These 
features can be attributed to the fact that detonation 
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waves with reaction zones governed by activated (Ar¬ 
rhenius) kinetics are unstable to perturbations in the 
transverse direction. m\ In condensed-phase explosives 
(i.e., polycrystalline high explosives), the reaction zone 
is usually controlled by the heterogeneity of the medium, 
wherein shock localization at density inhomogeneities 
results in local centers of energy release, so-called “hot 
spots,” from which reaction fronts burn-out the rest of 
the explosive within the reaction zone. Remarkably, 
despite these multidimensional and unsteady features, 
the average velocity of detonation fronts is usually very 
close to the predictions of the CJ criterion, provided 
that the detonation is not near its propagation limits. 
In gases, the average velocity of a detonation wave is 
typically within 1% of the ideal CJ velocity, despite the 
fact that the local shock front velocity may undergo 
transient deviations as large as 50% over the cycle of a 
detonation cell. 

The usual explanation offered to account for the 
success of the CJ criterion is that it can still be applied 
in an average sense to unsteady and multidimensional 
detonation waves. For example, if a large enough control 
volume, moving at a steady velocity, can enclose the 
detonation wave, then the CJ criterion (along with the 
conservation laws) can still be applied with accuracy to 
the average flow entering and exiting the control volume. 
This answer is not entirely satisfying, since the averaging 
procedure may not be justified. The motivation for the 
present paper is to put this CJ criterion to a rigorous 
test by considering an extreme case of a highly discrete 
media wherein all of the energy release is concentrated in 
(5-functions in space and time, embedded in an otherwise 
inert medium. This scenario can be considered the 
extreme limit of a heterogeneous energetic medium. 
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in contrast with the spatially uniform energetic media 
considered in the CJ model. 

A similar examination of reaction-diffusion waves in 
biological excitable media and flames in energetic media 
in the limit of spatially discrete sources has recently 
been undertaken. [HHni This approach has generated a 
number of unique results that would not be encountered 
in traditional flame theory, such as a limit to flame 
propagation in adiabatic systems that is distinct from 
the classical thermodynamic limit. [IMlj Experimental 
results showing an independence of flame speed on oxy¬ 
gen concentration for fuel particulates in suspension in 
a gaseous mixture has provided experimental validation 
of another prediction of discrete-source flame theory. [T5] 
This approach has resulted in recognition of a separate 
branch of combustion, i.e.,“discrete combustion,” that 
is unique from both traditional homogeneous and 
heterogeneous combustion. m The present paper is an 
initial examination into the analogous problem of a 
discrete-source detonation. 

Exploring the unsteady or multidimensional detona¬ 
tion dynamics usually necessitates numerical simulations 
due to the inherent nonlinearity of the governing Euler 
equations of compressible fluid flow. The ability to 
numerically resolve spatial and temporal d-functions 
is a challenging problem, and thus the examination of 
simpler, analog systems for which analytic solutions 
are available is of interest. The present paper utilizes 
a reactive Burgers equation analog that was originally 
formulated by Fickett [HHH] and Majda [H]. The 
advantage of using this analog is its simplicity. Instead 
of the conservation of mass, momentum, and energy in 
the Euler equations, the conservation of only one anal¬ 
ogous scalar quantity is considered in this analog. This 
simplification has been rigorously justified by Faria et 
al. showing that the reactive compressible Navier-Stokes 
equations can be reduced to a one-dimensional reactive 
Burgers equation based on an asymptotic analysis. [20] 

This analog approach to examining detonation 
dynamics has recently generated a variety of noteworthy 
results, including the ability to capture the period dou¬ 
bling sequence of bifurcations that leads to a chaotically 
pulsating detonation wave pDII23j . and to qualitatively 
model the phenomenon of shock-induced ignition in a 
spatially uniform reactive media governed by two-step 
induction-reaction kinetics [24] . This approach has been 
expanded to study non-ideal one-dimensional detona¬ 
tions and was able to capture the instabilities caused 
by the presence of a loss term that can represent the 
effects of front curvature or friction on the detonation 
wave.pS] All of these prior studies treated the medium 
as being spatially uniform. The present paper is the first 
to examine discrete-source detonations in analogs based 
upon the Burgers equation. 


II. DESCRIPTION OF ANALOG MODEL 
A. Burgers Analog Model 

Using the one-dimensional scalar Burgers equation as 
a model to explore detonation phenomena was proposed 
by Fickett [TMll^ and independently by Majda [12]. For 
this paper, we will construct our system in an Eulerian 
reference frame, following Fickett [TB], however, the en¬ 
ergy released by reaction will be included as a source 
term in the conserved variable, as in Majda’s model m- 
This analog system is based on the conservation of E, a 
variable representing some features of an extensive prop¬ 
erty, such as mass, momentum, or energy, in compressible 
flow. The volumetric density of A, denoted p, is analo¬ 
gous to an intensive property, such as density, velocity, or 
specific energy, in compressible flow. The governing con¬ 
servation equation is in an Eulerian reference frame, in 
other words, the conservation law is applied to a specific 
region of space through which fluid flows, i.e., a control 
volume. For a reactive flow, E is convected into and out 
of the control volume by the flow while also being gen¬ 
erated by chemical reaction within the control volume. 
The diffusion of E through the boundary of the control 
volume is neglected in this model, i.e., an inviscid flow is 
assumed. 


Control volume 
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FIG. 1. Conservation of E in the control volume 

In Fig. f{p) is the flux of E entering and leaving 
the control volume, which is analogous to pressure in the 
momentum equation. E^^ is the time rate of change of 
E in the control volume. The generation rate of energy 
E within the control volume is equal to the time rate of 
change of the reaction progress variable, Z, multiplied 
by Qcv! total energy released in the control volume 
when the reaction is complete. Note that the source 
energy is a property of space only, so it is not 
advected or compressed in this analog system. In other 
words, energy sources are assumed to be stationary in a 
lab-fixed reference frame. 

The equation of energy conservation for the control 
volume can be written as follows, 

4, =Q^,Z + /(p,)-/(p,+a.) (1) 

Dividing by Acc on both sides of Eq. [^ the equation 
of energy conservation can be obtained in terms of the 
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intensive property, p, 

• _ Qcv ^ / jPx) ~ / jPx+Ax) /2\ 

^ Act Ax 

Here, Z represents the mass fraction of product (i.e., 
Z = 0 for unreacted and Z = 1 for fully reacted) in this 
system. 

Introducing the spatial density of the energy source, 
q = , and taking the limit of Eq. as Ax approaches 

0, the following equation can be obtained, 

By introducing c = ^, Eq. can be written as 
ot ox 

Equation!^ is a nonlinear advection equation with reac¬ 
tive source term. Hence, c represents the speed at which 
p is advected. The solution p(x, t) is constant along each 
curve x—ct = Xg in the x-t plane. These curves are known 
as the characteristics of the equation, and c is the sound 
speed relative to a lab-fixed reference frame, which is the 
slope of the corresponding characteristic in the x-t plane. 


In general, a reaction law is required to specify 
the dependence of reaction rate, on Z and p. 

However, as the simplest approximation, the assumption 
of instantaneous reaction, i.e., an infinitely thin reaction 
zone, is made for the following analysis. Hence, no 
reaction law needs to be specified. 


B. Solution for Spatially Homogeneous Source 
Energy Release 

In a homogeneous reactive system, where the density 
of the energy source is uniform, g is a constant. By 
applying the CJ condition, that the fully reacted flow 
is exactly sonic relative to the leading shock, i.e., 
Tj = Pcj = the CJ detonation velocity for the 

homogeneous system with density of energy source q can 
be obtained as 

= 2g (9) 

The details of this derivation can be found in Ap¬ 
pendix 1^ 


The relation between / and p is considered as the 
equation of state. The function of f{p) must satisfy the 
entropy condition for a shock wave solution, which re¬ 
quires characteristics propagate into the shock in the x-t 
diagram, as time advances. [2(1) For a shock propagating 
rightwards, the shock speed S satisfies the entropy con¬ 
dition if 


dx 


Pl 


c^> S > 



( 5 ) 


where subscripts “R” and “L” indicate the pre-shock and 
post-shock states, respectively. Eq. requires f{p) to 
be a convex function, i.e., > 0, hence, a quadratic 

function is chosen as the simplest approximation of the 
equation of state in the detonation analog [T?>HTn] 



( 6 ) 


From Eq. a relation between c and p can be obtained. 


dp 


= c = p 


( 7 ) 


which means that the advection velocity of the local en¬ 
ergy density p equals to the value of p itself. Inserting 
Eq. into Eq. the one-dimensional inviscid reactive 
Burgers equation arises. 


di ^ ^ 


( 8 ) 


The solution of p for a homogeneous source energy 
is a right triangular profile over space with a height of 
p^j at any time. This triangle gets stretched over time, 
with its vertical leg propagating to the right at velocity 
The area under the triangular profile of p is the 
total amount of energy released over the distance traveled 
by the detonation wave at t. The solution of p for a 
homogeneous system is illustrated in Fig. 



FIG. 2. Sample solution for spatially homogeneous source 
energy release 


C. Solution for Discrete Source Energy Release 

For a heterogeneous system, the source energy density 
is not uniform over space. Hence, a spatial distribution 
of q, i.e., q{x), must be specified. In this study, the limit 
of a highly heterogeneous system is considered wherein 
the release is concentrated in point-like sources in space. 
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For this case, q{x) can be specified as a summation of 
equally or randomly spaced ^-functions, 

«(-) = E^(^) (10) 

2 ^ '' 

where A is the total amount of energy released by each 
discrete source, x^., is the location of each discrete 
source, and A^totai is the total number of energy sources 
within an arbitrary length of the reactive system, Ftotai- 
Hence, given A, A^totai, and Ltotai, the average energy 
density of a discrete system, (^discrete, can be determined. 
The average energy density is equivalent to the uniform 
energy density of a homogeneous system. 


Q^discrete 



Etotal 


ANi 


total 


A 


total 


Etotal -^total 

Recall that S{x) dx = 1. It is important to note that, 
for a given A, N-^, and Ltotai, the discrete systems with 
both equally and randomly spaced energy sources have 
the same global energy release as the homogeneous sys¬ 
tem. However, for a discrete system with equally spaced 
energy sources, this equivalence can be expressed in a 
simpler form 


^discrete — ^ — Q 


( 12 ) 


where L is the spacing between each two adjacent energy 
sources. The source can have a finite delay period r 
between being shocked and releasing energy. 


Except for the instants in which the sources are acti¬ 
vated, the solution of the discrete system is governed by 
the non-reactive Burgers equation. 


dp dp 
dt^^dx 


= 0 


(13) 


As this system has only right-running characteristics, the 
solution of a shock propagating with a linear profile of p 
upstream and downstream of the wave can be treated 
analytically using the method of characteristics. The de¬ 
tailed solution for the propagating shock can be found in 
Appendixj^ The shock-shock and shock-contact interac¬ 
tions can be similarly treated analytically. A schematic 
x-t diagram of how the method of characteristics is used 
to solve for the trajectory of the shock waves, contact dis¬ 
continuities, and their subsequent interactions is shown 
in Fig.d 


III. IMPLEMENTATION OF ANALOG MODEL 

The methodology of constructing a solution to the 
problem of a shock propagating in a system consisting of 



FIG. 3. Schematic illustration of the method of characteristics 
in x-t diagram. The right-running characteristics are plotted 
as dashed, arrowed lines, the trajectories of contact discon¬ 
tinuities as solid, arrowed lines, the trajectories of shocks as 
thick solid curves, and the locus of sources releasing energy as 
black solid circles that become red open squares upon being 
shocked and red solid squares upon releasing their energy. 


discrete sources is implemented via a computer program 
that tracks the trajectory of the leading shock, internal 
shocks, contact discontinuities, and the slopes of the saw¬ 
tooth regions between these various elements. Although 
a computer is used to track the interacting elements, we 
emphasize that the solution is entirely analytic and no 
finite difference approximation is made. When the lead¬ 
ing shock has crossed a source location and the prescribed 
delay period has elapsed, the solution is paused and the 
new source inserted into the profile of p at the source lo¬ 
cation. To implement this, the d-function, with an area 
of A under it, is approximated as a triangular profile with 
very small width w and height h =—, as illustrated in 
Fig. a Sample solutions of p for randomly discrete sys¬ 
tems are shown in Fig. 

The only convergence parameters in this method are 
the source width w and the precision with which the time 
of element interactions (shock-shock, shock-contact) are 
solved for via a numerical root-finding algorithm. It is 
relatively easy to verify that the solution is independent 
of these parameters once they are made sufficiently small 
(width of ^-function w/L = 10“^° and convergence fac¬ 
tor in solving for time of wave intersections eltc = 10“^, 
where 4 = L/D^^ is the characteristic time for this 
study). This convergence study was performed in the 
range of w/L = 10“®-10“^^ and e/t^ = 10“^-10“® for 
all results reported in this paper, and the results were 
found to be independent of the value of these parameters. 


IV. RESULTS 


In the results reported below, the spatial coordinate 
X is normalized by the average spacing between each two 
adjacent sources L, leading shock velocity D by and 
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FIG. 4. Illustration of the implementation of source activa¬ 
tion: (a) source during delay period after being shocked, (b) 
source energy release inserted as a triangular profile to ap¬ 
proximate (5-function, and (c) shortly after being activated, 
and (d) the second shock now caught up and merged with 
the leading shock. 


time t by tc- Thus, the dimensionless results are only 
dependent upon q and r. 


A. Equally Spaced Energy Sources with Zero Delay 


For the case of zero delay, i.e., a source releases its 
energy immediately upon contact with the leading shock 
front, it is possible to solve analytically for the average 
propagation velocity without having to rely upon a com¬ 
puter program to track the wave interactions. In this 
case, the ^-function is inserted at the front and there is 
no interaction between the resulting blast wave and the 
pre-existing sawtooth profile. In this case, the decay of 
the shock front is given by m 


= \p2At 


(14) 


From Eq. 14 the time required for the shock front to 


reach the next source a distance L away can be solved 



(15) 



X 


FIG. 5. Sample solutions of p for randomly spaced energy 
sources with random period delay. Circles are unshocked 
sources, empty squares shocked sources during delay period, 
and solid squares sources with energy released. 


Hence, the average velocity for the shock front to travel 
to the next source over a distance L is 


_ L 

^avg — T| 


2A 

~L~ 


2g = -Dcj 


(16) 


which is equal to the CJ detonation velocity of the equiv¬ 
alent homogeneous media. The detailed derivation of 
= 2q can be found in Appendix Since there is 
no influence from one cycle to the next, the wave will 
continue to propagate in a periodic manner, with the av¬ 
erage wave speed always equal to this value. 


B. Equally Spaced Energy Sources with Fixed 
Period Delay 

The results shown in Fig. are of a case where 
discrete sources were equally spaced along the x-axis 
with an average energy density of g = 0.5 (i.e., = 1) 

and fixed delays of r = 0.5tc and r = O.ltc- 

In Fig.ga), the history of the instantaneous and av¬ 
erage leading shock velocities for r = 0.5tc are compared 
to the CJ velocity of the system. The instantaneous 
velocity was obtained analytically from the solution, 
while the average velocity was computed by simply 
dividing the distance the leading shock had traveled by 
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the time elapsed. 




xlL xlL 


FIG. 6. Velocity history of a system where sources were 
equally spaced with an average energy density of q — 0.5 and 
(a) a fixed delay of r = 0.5tc over a short distance in linear 
scale, and (b) fixed delays of r = O.SG and r = O.IG over a 
long distance in logarithmic scale. 


As shown in Fig. [^a), the instantaneous velocity 
fluctuated between values higher and lower than the 
CJ velocity, while the average velocity slowly ap¬ 
proached the CJ velocity. This fluctuation appeared 
to be in self-repetitive cycles as the average velocity 
approached the CJ value. As shown Fig. I^b), for a 
short delay time, t = O.ltc, the wave is overdriven early 
in the propagation. As the leading shock traveled a 
sufficiently long distance {x/L = 200), the average veloc¬ 
ities for both cases of r = 0.5tc and r = O.ltc converge 
to the CJ velocity with a difference less than 2.5% of . 

To better examine the interaction of waves and source 
energy release during the propagation of a detonation 
wave, the location and time of source activation, and the 
trajectories of shocks and contact surfaces were plotted 
in the x-t diagrams shown in Fig. [^a) and (b). 

Figure [^a) shows the x-t diagram constructed in a 
lab-fixed reference frame. A shock and a contact surface 
originate each time a new source releases its energy. The 
characteristics are not shown in this figure for clarity, but 
they have a linear profile between the shocks and contact 
surfaces, as shown in the schematic of Fig. The shock 
generated by a newly released source eventually catches 
up and accelerates the leading shock. It is interesting to 
note that, as the detonation propagates, the trajectories 
of the contact surfaces, i.e., the trailing characteristics of 
the rarefaction waves, tend to coalesce along the path in 
x-t space where source energy release events are occur¬ 
ring. 

Figure Hb) shows the x-t diagram constructed in a 
reference frame moving at the CJ velocity of a detonation 
in the equivalent homogeneous system. The time interval 
{t = 40 to 45tc) shown in Fig. [^b) is after the leading 
shock has traveled a sufficiently long distance so the av¬ 
erage velocity is very close to the CJ value. The leading 
shock front periodically oscillates around approximately 



(x-D^,t)/L 


(x-D^,l)IL 


FIG. 7. x-t diagrams of wave processes, (a) in a lab-fixed ref¬ 
erence frame and (b) in a wave-attached frame moving at the 
GJ velocity, of a system where sources were equally spaced 
with an average energy density of g = 0.5 and a fixed delay of 
T = 0.5tc. x-t diagrams of wave processes, (c) in a lab-fixed 
reference frame and (d) in a wave-attached frame moving at 
the GJ velocity, of a system where sources were randomly 
spaced with an average energy density of q = 0.5 and random 
delay in the range 0.1 to It^. The trajectories of contact dis¬ 
continuities are plotted as blue lines, the trajectories of shocks 
as black curves, and the locus of sources releasing energy as 
black solid circles that become red open squares upon being 
shocked and red solid squares upon releasing their energy. 


the same location relative to the reference frame moving 
at . The sources release their energy at almost equal 
time intervals. As can be clearly observed in this figure, 
the trailing characteristics associated with prior sources 
pass extremely close to the locus of new sources. 


C. Randomly Spaced Energy Sources with 
Random Period Delay 

The locations of energy sources are distributed within 
a distance Ttotai using a random number generator, i.e., 
~ U[0, Ltotai], and the delay time associated with 
each source is randomly generated within a prescribed 
range, i.e., ^ U[r„iin, Tmax]- The results shown in this 

section are for the case where discrete sources of energy 
were randomly spaced with an average energy density of 
q — 0.5. In this case, L is the average spacing between 
sources, L = . 
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FIG. 8. (a) Velocity history of a system where sources were randomly spaced with an average energy density of q = 0.5 and 
random time delays over the range 0.1 to Itc over a short distance in linear scale, (b) The average wave velocity for a set of 
20 realizations of a system of randomly spaced sources with an average energy density of q = 0.5 and random time delays over 
the range 0.1 to Itc over a long distance in logarithmic scale, (c) The average wave velocity for two sets of 20 realizations of 
systems of randomly spaced sources with an average energy density of q = 0.5 and random time delays over the ranges 0.01 to 
O.ltc and 1 to lOG over a long distance in logarithmic scale. The thick dashed lines in (b) and (c) are the ensemble averaged 
velocities as a function of x/L. 


Figure [^b) shows the average wave velocity as a 
function of propagation distance for a set of 20 real¬ 
izations of a system with randomly spaced sources and 
random time delays over the range 0.1 to Itc- While the 
average velocities exhibit wide fluctuations, particularly 
early in the propagation, the average velocities all 
approach the equivalent CJ velocity of the homogenized 
system. The fact that many of the individual realizations 
exhibit very large velocities early on is attributed to the 
fact that random placement of sources makes it probable 
that two sources will be located quite close together, 
resulting in an initially highly overdriven wave. A thick 
dashed line, corresponding to an ensemble average of 
20 such solutions, is also shown to assist the reader 
in identifying the overall trend. In Fig. [^c), similar 
ensembles are shown for the cases of random delay 
times over the interval 0.01 to O.ltc and the interval 
1 to lOtc- In the case of the very short delay, the 
ensemble averaged velocity decays monotonically to 
the CJ velocity of the homogenized system. For the 
case with long delay time, most of the individual cases 
decay to well below the equivalent CJ velocity due to 
the time required for the shock waves generated by 
the sources to catch up to the leading shock. Upon 
encountering hundreds of subsequence sources, however, 
the average velocity exhibits a steadily increasing value, 
asymptotically approaching the CJ velocity. In all cases, 
after propagating through 10^ sources, the average wave 
velocity is within 1% of the CJ velocity. 

The wave processes of this system with randomly 
spaced sources and random delay are plotted in the x-t 
diagram shown in Fig. [^c) and (d). The x-t diagrams 
constructed in both a lab-fixed reference frame and 
a reference frame moving at are shown. Again, 
the characteristics in the regions of constant slope 


between the various shocks and contact surfaces are not 
plotted in these figures for clarity. No self-repetitive 
wave interactions can be observed even after the wave 
front had traveled a significantly long distance. The 
phenomenon of the trailing characteristics of rarefaction 
waves coinciding the loci of source firings that was ob¬ 
served with a regular spacing of sources (Section IV B) 
is not observed in the case with random sources. 


V. DISCUSSION 

When the release of the conserved variable (here, rep¬ 
resenting energy) is concentrated in discrete, point-like 
sources and then released when initiated by a passing 
shock wave in the system described by the Burgers 
equation considered in this paper, the resulting blast 
waves interact and merge with each other and with the 
leading shock. While the instantaneous velocity of this 
front exhibits large variations, the average velocity of 
the front after traversing many hundreds or thousands 
of sources is very close (within 2.5%) of the ideal 
Chapman Jouguet velocity of the equivalent uniform 
medium in which the sources have been homogenized 
throughout the medium. This result was observed both 
with regularly spaced sources with a fixed delay time 
between the shock and with systems with randomly 
distributed sources and random delay time, as can 
be seen in Figs. and respectively. This result 
provides a strong validation of the CJ criterion for 
highly spatially heterogeneous reactive media, at least 
within the Burgers equation analog system. It may also 
provide some rationale as to why the CJ criterion is so 
successful in describing detonations in gases which are 
characterized by an unstable, multidimensional cellular 











structure. If this result, namely that a detonation will 
always propagate on average at the velocity of the CJ 
detonation in an equivalent homogeneous system, is 
proven to be universal, it may even supplant the tradi¬ 
tional CJ criterion, which was formulated for detonation 
waves propagating at steady velocity in uniform media. 

When a source of the conserved variable is released as 
a J-function, it evolves into a triangular, sawtooth profile 
with a leading shock followed by a linear rarefaction. In 
the present study, the source release was initialized as 
a triangular profile, but it can be shown that any non¬ 
negative profile with a local maximum (i.e., single hump) 
will eventually develop into this forward propagating 
triangular (or sawtooth) profile wave |27) . Thus, we do 
not believe our results are influenced by the selection of a 
triangular profile in introducing the sources. The arrival 
of the shock at the leading front, and the subsequent 
merging of these shocks, accelerates and sustains the 
leading front, while the following rarefaction results in 
a continued deceleration of the wave. For the case of 
regularly spaced sources with a fixed delay period, the 
time at which the sources release their energy evolves 
to become very close to the time of arrival of the locus 
of trailing rarefactions from all prior sources, as seen in 
Fig.0 This phenomenon, in which the trailing rarefac¬ 
tions appear to pile-up on top of each other and become 
coincident with the time at which the new sources are 
released, resulted in considerable numerical difficulty in 
implementing the event-driven solution technique used 
in this paper, since it was necessary to determine the 
intersection point in space and time of all interacting 
elements (shocks, trailing rarefactions, and new sources) 
in order to ascertain which interaction occurs next. This 
condition may also be the discrete-source analog to the 
traditional CJ criterion, wherein the leading edge of 
the trailing rarefaction of the unsteady expansion of 
detonation products becomes coincident with the sonic 
locus at the end of the reaction zone. 

The treatment of detonation wave dynamics by 
considering a detonation as an ensemble of interacting 
waves originating from discrete sources may have ap¬ 
plication beyond just computing the average velocity of 
detonation propagation in an effectively infinite medium. 
When losses are present (e.g., lateral expansion due to 
yielding confinement, momentum losses due to boundary 
layers in narrow channels, etc.), detonations no longer 
propagate at the ideal CJ detonation velocity even 
in homogeneous media, and calculation of the correct 
propagation velocity necessitates solving for the reaction 
zone structure with losses present. The calculation 
procedure consists of solving an eigenvalue problem 
with the post-shock state forming one boundary and a 
reaction zone structure that passes smoothly through 
a sonic point (i.e., saddle condition) defining the other 
boundary. While this technique is well defined for a 
smooth, laminar-like wave structure in a homogeneous 


medium (see Ref. [IH]), the procedure for heterogeneous 
or multiphase detonations has yet to be formulated, and 
must be treated via direct numerical simulation. Results 
in recent years examining detonations in condensed- 
phase explosives with large-scale heterogeneities have 
generated results that cannot be explained using 
conventional front curvature theory, which assumes a 
laminar-like wave structure. In the presence of losses, 
a system of interacting, multidimensional blast waves 
originating from spatially discrete sources may give 
rise to a percolating-like phenomenon in which waves 
would be able to propagate from source to source 
under conditions wherein a wave in the equivalent 
homogeneous medium would fail. This type of be¬ 
havior has been recently demonstrated in detonation 
simulations using the reactive Euler equations with a 
medium in which a spatial inhomogeneity has been 
introduced. [29] This approach to detonation dynamics 
may provide a theoretical framework to account for the 
anomalous results in highly heterogeneous explosives. |5D] 

Since the resulting wave structure in the system stud¬ 
ied in this paper is seen to be an interacting ensemble of 
decaying sawtooth-profile blast waves described by the 
classical inviscid Burgers equation, it would be of inter¬ 
est to further explore these results within the framework 
of Burgers turbulence. In recent years, there has been 
a resurgence of examination of the Burgers equation to 
study scaling exponents in the power law of the proba¬ 
bility distribution as solutions of the equation, initialized 
with white noise for example, decay to a sawtooth profile 
(see thorough review by Bee and Khanin jUj). The fact 
that turbulence-like scaling is obtained from the inviscid, 
scalar Burgers equation, wherein all viscous dissipation 
occurs in the discontinuities embedded within the solu¬ 
tion, is a remarkable finding. Some of the results of the 
present study, in which randomly spaced sources were 
triggered with random delay times, could be considered 
in the class of kicked Burgers turbulence, but under a 
scenario wherein the forcing function is triggered or syn¬ 
chronized with the passage of the shock front. 


VI. CONCLUSION 

Solutions to the inviscid Burgers equation with spa¬ 
tially and temporally discrete sources of the conserved 
variable are released when activated by passage of the 
leading shock front were examined. The governing equa¬ 
tion was non-reactive outside of the instants in time when 
the release of the sources occurred, enabling the solution 
to be constructed entirely analytically. The release of 
a source was treated by inserting a J-function into the 
solution, and then allowing the resulting waves to inter¬ 
act. The subsequent interaction of the blast waves in the 
Burgers equation was observed to propagate, on average, 
at a speed identical to the Chapman-Jouguet speed of a 
steady detonation in an equivalent homogeneous system 
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with the same average spatial release of the conserved 
variable. This result applied to both systems with uni¬ 
formly spaced sources and fixed periods of delay between 
shock and release as well as systems with randomized 
energy release (i.e., sources randomly spaced and with 
random delay time, but with the sources still activated 
by the passage of a shock). In the regular spacing and 
delay case, as the solution progressed, the locus of trail¬ 
ing rarefaction waves evolved to become coincident with 
the release of new sources, a result that may have some 
connection to the classical CJ criterion for homogenous 
media. In a randomly spaced system with random de¬ 
lays, the release of the sources and rarefactions from prior 
sources did not exhibit this coincident behavior, but the 
resulting wave complex still propagated with an average 
velocity equal to that predicted by the CJ criterion of the 
equivalent homogeneous medium. 

Appendix A: Chapman-Jouguet Detonation Velocity 
for a System with Spatially Homogeneous Energy 
Release 

For a control volume enclosing and attached to a non¬ 
reactive shock wave (as illustrated in Fig. [^, the flux of 
E leaving the control volume must equal to that entering 
the control volume 


This is the well-known result that a shock, as a solution 
of the Burgers equation, propagates at the average of 
the pre-shock and post-shock characteristic speed. 

For a detonation wave, the difference between the 
fluxes leaving and entering the control volume equals the 
rate of E being released or absorbed by reaction within 
the control volume. Since the energy source is assumed 
to be stationary in a lab-fixed reference frame, it enters 
and leaves the control volume at the speed of detonation 
propagation, D. Thus, the conservation equation across 
a detonation wave can be written as 

{D - U 2 )p 2 - {D - ui)pi = D {Z 2 - Zi)q (A4) 

Note that Z\ is always equal to 0 since the reactant 
in front of the detonation wave is fully unreacted. By 
/(p) = up, D can be expressed as 

D = (A5) 

(P2 - qZ2) - Pi 

For any value of Z 2 such that 0 < Z 2 < 1, Eq. |A5| can be 
plotted as a straight line with slope D in p-f{p) space, 
which is analogous to the Rayleigh line for the governing 
equations of compressible fluid flow. The analog of the 
Hugoniot curve is Eq. the equation of state. 


(5 - U 2 ) P2- {S - ui) Pi = 0 


(AI) 


where S is the speed of shock propagation, and u is the 
flow speed with respect to a lab-fixed reference frame. 
Since /(p) = up, S can be derived from Eq. AI as 


The Chapman-Jouguet condition, which states that 
the flow becomes sonic relative to the detonation, i.e., 
Tj = Pcj ~ when the material is fully reacted, i.e., 
^2 = I, requires the following relations to be satisfied 
simultaneously. 


/(P 2 ) - f{pi) 
P2 - Pi 


(A2) 


which is known as the Rankine-Hugoniot jump condition. 


Control volume 




Here, pi = 0 as the leading shock wave is assumed to 
be very strong, and the equation of state has the form 
of Eq. By solving Eq. |A6[ the CJ detonation speed 
for a homogeneous system with energy density q can be 
obtained as 


Da, = 2g (A7) 

As shown in Fig. ^ can be equivalently de¬ 
termined by finding the slope D corresponding to the 
Rayleigh line with Z 2 = 1 tangent to the Hugoniot curve 
f{p) = \P^ in p-/(p) space. 


FIG. 9. Schematic representation of the control volume en¬ 
closing a nonreactive shock wave. 

Considering the specific equation of state as Eq. S 
can be solved as 

S =]^{pi+ P 2 ) (AS) 


Appendix B: Propagation of a Shock Wave Solved 
via the Method of Characteristics 

In the system with energy sources consisting 5- 
functions in space, the shock waves always propagate 
with a linear profile of p upstream and downstream of 
the shocks. For a special case of the leading shock with a 












(a) 
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FIG. 10. The analog of the Chapman-Jouguet condition 
shown in p-f{p) diagram. 


zero-slope profile of p upstream, the position of the shock 
front as a function time can be solved as Eq. Fig¬ 

ure [Ufa) is a snapshot of a shock located between two 
linear profiles of p at an arbitrary time. For simplicity, 
we set the location of the shock to be x = 0 at this instant 
t = 0. The values of p immediately in front and behind 
the shock are pi and p 2 , respectively. The upstream and 
downstream linear profiles of p at t = 0, i.e., po,i and po ,2 
can be expressed as follows. 

Pop = Pi + Po ,2 = P 2 + X 0 S 2 (Bl) 

where si and S 2 are the slopes of the upstream and 
downstream linear profiles, respectively, and xq is the 
location of each characteristic at t = 0. 

In the proposed system, since signals propagate at a 
sound speed equal to the corresponding value of p and 
only in the rightward direction, the characteristics in an 
x-t diagram are straight lines along which p is constant 
and always equals to its initial value at t = 0. The func¬ 
tion of each characteristic in upstream and downstream 
regions is as follows. 


Xi = Xo -I- popt X2= xq+ po,2t (B2) 

Combining Eqs. [BT]and[B^ the constant value of p along 
each characteristic can be expressed as 


Po,i — 


SlXi + Pi 

S\t -\- 1 


P 0.2 — 


S2X2 P2 
S2t 1 


(B3) 


As shown in Fig. the shock propagates along a 

trajectory where the characteristics from the upstream 
region intersect with those from the downstream region 
in x-t space, determining the shock velocity. As obtained 
in Eq. |A3[ the speed of shock propagation at any time is 
equal to the average of the sound speeds corresponding 
to the two characteristics intersecting at this time. Since 
along each characteristic, the sound speed always equals 
to the initial value of p at t = 0, Eq. |A3|can be rewritten 


@ ® 


Poi= P2 + X0S2 .P^ 





(b) 



FIG. 11. Illustration of (a) a shock wave located between two 
linear profiles of p upstream and downstream of the shock and 
(b) the trajectory of the shock wave solved using the method 
of characteristics in an x-t diagram. 


as 




(B4) 


Using Eq. B3 to express po,i and po ,2 in Eq. B4 and 
substituting xi and X 2 by the position of the shock, Xg, a 
first-order linear ordinary differential equation governing 
the propagation of the shock wave is obtained. 


S = 


dt 


1 f siXi -I- Pi 522:2 + P2 


2 y S\t -\- 1 


S2t -f 1 


(B5) 


with initial condition x^ = 0 at t = 0. By integrating 
Eq. |B5| the position of the shock as a function of time 
can be solved, 

2 ^S (0 = s Is [pi (VSit + l-\/ S 2 t -\- 1 — S 2 t — 1 ) 


- P 2 (VsU -b lVs 2 t -b 1 - Sit - 1 ) ] (B6) 

This solution is used in the main paper to track the tra¬ 
jectory of the shock waves. 
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